;;; Exploration and algorithms for fractals and fractal landscape generation

(in-package :cl-user)

;;; Mid Point Displacement Algorithm

(defun basic-displacement-function (&optional (max-step 2))
  (if (= (random 2) 1)
      (random max-step)
      (- (random max-step))))

(defun mpd (a b c d steps &optional (displacement-function #'basic-displacement-function))
  "Create a height array using the mid point displacement algorithm.  Initial heights are in a, b, c, d, and the algorithm runs number of times specified in steps."
  (let* ((max-xy (1+ (expt 2 steps)))
	 (dims (list max-xy max-xy))
	 (heights (make-array dims)))
    (flet ((mpd-step (a b c d)
	     (let* ((e (+ (/ (+ a b c d) 4) (funcall displacement-function)))
		    (f (+ (/ (+ a c e e) 4) (funcall displacement-function)))
		    (g (+ (/ (+ a b e e) 4) (funcall displacement-function)))
		    (h (+ (/ (+ b d e e) 4) (funcall displacement-function)))
		    (i (+ (/ (+ c d e e) 4 ) (funcall displacement-function)))))))
      (setf (aref heights 0 0) a
	    (aref heights max-xy 0) b
	    (aref heights 0 max-xy) c
	    (aref heights max-xy max-xy) d)
      heights)))


;;;; Diamond-Square Algorithm (http://www.gameprogrammer.com/fractal.html#other)

(defun random-range (min max)
  (let ((scale (- max min)))
    (+ (random scale) min)))

;;; Given centerpoint C, returns the values in the height map
;;; as indicated by v in the diagram below.
;;; . v .
;;; v C v
;;; . v .
(defun diamond-adjacent-grid-values (heights i j i_len j_len)
  (let ((imax (1- (array-dimension heights 0)))
	(jmax (1- (array-dimension heights 1)))
	(i+ (+ i i_len))
	(i- (- i i_len))
	(j+ (+ j j_len))
	(j- (- j j_len)))
    (if (< i- 0) ; top row
	(if (< j- 0) ; Top left
	    (values (list (aref heights i+ j)
			  (aref heights i j+))
		    2)
	    (if (> j+ jmax) ; top right
		(values (list (aref heights i+ j)
			      (aref heights i j-))
			2)
		(values (list (aref heights i+ j)
			      (aref heights i j+)
			      (aref heights i j-))
			3)))
	(if (> i+ imax) ; bottom row
	    (if (< j- 0) ; bottom left
		(values (list (aref heights i- j)
			      (aref heights i j+))
			2)
		(if (> j+ jmax) ;; bottom right
		    (values (list (aref heights i- j)
				  (aref heights i j-))
			    2)
		    (values (list (aref heights i- j)
				  (aref heights i j-)
				  (aref heights i j+))
			    3)))
	    (if (< j- 0) ; left edge
		(values (list (aref heights i- j)
			      (aref heights i+ j)
			      (aref heights i j+))
			3)
		(if (> j+ jmax) ; right edge
		    (values (list (aref heights i- j)
				  (aref heights i+ j)
				  (aref heights i j-))
			    3)
		    (values (list (aref heights i- j) ;; Somewhere in the middle
				  (aref heights i+ j)
				  (aref heights i j-)
				  (aref heights i j+))
			    4)))))))

(defun test-array ()
  (make-array '(3 3) :initial-contents '((1 2 3) (4 5 6) (7 8 9))))

(defun bigger-test-array ()
  (let ((arr (make-array '(9 9))))
    (dotimes (i 9 arr)
      (dotimes (j 9)
	(setf (aref arr i j) (+ j (* i 9)))))))

(defun test-dagv-1 ()
  (assert (null (set-difference '(2 4 6 8)
				(diamond-adjacent-grid-values (test-array) 1 1 1 1)))))

(defun test-dagv-corners ()
  (assert (null (set-difference '(2 4) (diamond-adjacent-grid-values (test-array) 0 0 1 1))))
  (assert (null (set-difference '(6 2) (diamond-adjacent-grid-values (test-array) 0 2 1 1))))
  (assert (null (set-difference '(4 8) (diamond-adjacent-grid-values (test-array) 2 0 1 1))))
  (assert (null (set-difference '(6 8) (diamond-adjacent-grid-values (test-array) 2 2 1 1)))))

(defun test-dagv-edges ()
  (assert (null (set-difference '(1 3 5) (diamond-adjacent-grid-values (test-array) 0 1 1 1))))
  (assert (null (set-difference '(1 5 7) (diamond-adjacent-grid-values (test-array) 1 0 1 1))))
  (assert (null (set-difference '(7 5 9) (diamond-adjacent-grid-values (test-array) 2 1 1 1))))
  (assert (null (set-difference '(3 5 9) (diamond-adjacent-grid-values (test-array) 1 2 1 1)))))

(defun test-dagv-3 ()
  (let ((arr (bigger-test-array)))
    (diamond-adjacent-grid-values arr 4 4 4 4)))

(defun test-dagv-4 ()
  (diamond-adjacent-grid-values (test-array) 0 0 1 1))


;;; Given centerpoint C, returns the values in the height map
;;; as indicated by v in the diagram below.
;;; v . v
;;; . C .
;;; v . v
(defun square-adjacent-grid-values (heights i j i_len j_len)
  (let ((imax (array-dimension heights 0))
	(jmax (array-dimension heights 1))
	(i+ (+ i i_len))
	(i- (- i i_len))
	(j+ (+ j j_len))
	(j- (- j j_len)))
    (if (< i- 0) ;; top row
	(if (> j+ jmax) ;; top right
	    (values (list (aref heights i+ j-)) 1)
	    (if (< j- 0) ;; Top left
		(values (list (aref heights i+ j+)) 1)
		(values (list (aref heights i+ j-)
			      (aref heights i+ j+))
			2)))
	(if (> i+ imax) ;; bottom row
	    (if (> j+ jmax) ;; Bottom right
		(values (list (aref heights i- j-)) 1)
		(if (< j- 0) ;; bottom left
		    (values (list (aref heights i- j+)) 1)
		    (values (list (aref heights i- j-)
				  (aref heights i- j+))
			    2)))
	    (if (> j+ jmax) ;; right edge
		(values (list (aref heights i- j+)
			      (aref heights i+ j+)
			      2))
		(if (< j- 0) ;; left edge
		    (values (list (aref heights i- j-)
				  (aref heights i+ j-)
				  2))
		    (values (list (aref heights i+ j-) ;; Somewhere in the middle
				  (aref heights i+ j+)
				  (aref heights i- j-)
				  (aref heights i- j+))
			    2)))))))
	
	

(defun do-step (heights i0 j0 i1 j1 displace)
  (let* ((i_len (truncate (/ (- i1 i0) 2)))
	 (j_len (truncate (/ (- j1 j0) 2)))
	 (ic i_len)
	 (jc j_len))
    ;; Diamond step
    (multiple-value-bind (vals cnt) (square-adjacent-grid-values heights ic jc i_len j_len)
      (setf (aref heights ic jc) (+ (/ (apply #'+ vals) cnt) (random-range (- displace) displace))))
    ;; Square step
    (multiple-value-bind (vals cnt) (diamond-adjacent-grid-values heights (- ic i_len) jc i_len j_len)
      (setf (aref heights (- ic i_len) jc) (+ (/ (apply #'+ vals) cnt) (random-range (- displace) displace))))
    (multiple-value-bind (vals cnt) (diamond-adjacent-grid-values heights (+ ic i_len) jc i_len j_len)
      (setf (aref heights (+ ic i_len) jc) (+ (/ (apply #'+ vals) cnt) (random-range (- displace) displace))))
    (multiple-value-bind (vals cnt) (diamond-adjacent-grid-values heights ic (+ jc j_len) i_len j_len)
      (setf (aref heights ic (+ jc j_len)) (+ (/ (apply #'+ vals) cnt) (random-range (- displace) displace))))
    (multiple-value-bind (vals cnt) (diamond-adjacent-grid-values heights ic (- jc j_len) i_len j_len)
      (setf (aref heights ic (- jc j_len)) (+ (/ (apply #'+ vals) cnt) (random-range (- displace) displace))))))

(defun random-array (rank &optional (max-random 1.0))
  (let ((a (make-array (list rank rank))))
    (dotimes (i rank a)
      (dotimes (j rank)
	(setf (aref a i j) (random max-random))))))

(defun ds-recurse (map i0 j0 i1 j1 displace)
  (let* ((i_side (truncate (/ (- i1 i0) 2)))
	 (j_side (truncate (/ (- j1 j0) 2)))
	 (i_c (+ i0 i_side))
	 (j_c (+ j0 j_side)))
    (unless (or (< i_side 1)
		(< j_side 1))
      (do-step map i0 j0 i1 j1 displace)
      (ds-recurse map i0 j0 i_c j_c displace)
      (ds-recurse map i_c j0 i1 j_c displace)
      (ds-recurse map i0 j_c i_c j1 displace)
      (ds-recurse map i_c j_c i1 j1 displace))))
  
(defun diamond-square-heightmap (map-power &key (max-displacement 1.0) (iterations 10) (roughness 1.0))
  (let* ((span (1+ (expt 2 map-power)))
	 (a (random-array span))
	 (disp max-displacement))
    (dotimes (n iterations a)
      (ds-recurse a 0 0 span span disp)
      (setf disp (* disp (expt 2 (- roughness)))))))
    
(defun pp-array (a)
  (dotimes (ri (first (array-dimensions a)))
    (format t "[ ")
    (dotimes (rj (second (array-dimensions a)))
      (format t "~A " (aref a ri rj)))
    (format t "]~%")))

(defun test-ds-recurse-1 ()
  (let ((map (random-array 5)))
    (pp-array map)
    (ds-recurse map 0 0 5 5 1.0)
    (pp-array map)))

;;;; Fluid Simulation (http://www.gamedev.net/reference/articles/article2001.asp)



